<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Example 6.9: Bounding consumer preference</title>
<link rel="canonical" href="/Users/mcgrant/Repos/CVX/examples/cvxbook/Ch06_approx_fitting/html/preference_regions.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Example 6.9: Bounding consumer preference</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 6.5.5, Figures 6.25-6.26</span>
<span class="comment">% Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX Argyris Zymnis - 11/30/2005</span>
<span class="comment">%</span>
<span class="comment">% We are given a set of consumer preference data for bundles</span>
<span class="comment">% of two goods x_1 and x_2. These points are generated by</span>
<span class="comment">% taking 40 random points and comparing them using the</span>
<span class="comment">% utility function: u(x_1,x_2) = (1.1*sqrt(x_1)+0.8*sqrt(x_2))/1.9</span>
<span class="comment">% Then, if we have u(i) &gt;= u(j) we say that (i,j) is in Pweak.</span>
<span class="comment">%</span>
<span class="comment">% Given this, we wish to compare the point (0.5,0.5) to each</span>
<span class="comment">% of the bundles in the given dataset. I.e. for each point k in the</span>
<span class="comment">% dataset, we wish to decide wether u(k) &gt;= u(0) or u(k) &lt;= u(0),</span>
<span class="comment">% or both, in which case we cannot make any conclusions about</span>
<span class="comment">% consumer preferences.</span>
<span class="comment">%</span>
<span class="comment">% To do this, we have to solve two LPs for each point:</span>
<span class="comment">%       minimize     u(k) - u(0)</span>
<span class="comment">%       subject to   g_i  &gt;= 0</span>
<span class="comment">%                    u(j) &lt;= u(i) + g_i^T(a_j - a_i), for all i,j</span>
<span class="comment">%                    u(i) &gt;= u(j), for all (i,j) in Pweak</span>
<span class="comment">%</span>
<span class="comment">% and:</span>
<span class="comment">%       maximize     u(k) - u(0)</span>
<span class="comment">%       subject to   g_i  &gt;= 0</span>
<span class="comment">%                    u(j) &lt;= u(i) + g_i^T(a_j - a_i), for all i,j</span>
<span class="comment">%                    u(i) &gt;= u(j), for all (i,j) in Pweak</span>
<span class="comment">%</span>
<span class="comment">% If the second LP has a strictly negative solution, we can deduce that</span>
<span class="comment">% u(k) &lt; u(0). If on the other hand the second LP has a nonnegative</span>
<span class="comment">% solution and the first LP has a strictly positive solution, we can</span>
<span class="comment">% deduce that u(k) &gt; u(0). Finally if none of the two previous cases</span>
<span class="comment">% holds, we cannot make a decision between the two bundles.</span>
<span class="comment">%</span>
<span class="comment">% NOTE: This file requires the auxilliary function utilfun.m to run.</span>

data= [<span class="keyword">...</span>
   4.5e-01   9.6e-01
   2.1e-01   3.4e-01
   9.6e-01   3.0e-02
   8.0e-02   9.2e-01
   2.0e-02   2.2e-01
   0.0e+00   3.9e-01
   2.6e-01   6.4e-01
   3.5e-01   9.7e-01
   9.1e-01   7.8e-01
   1.2e-01   1.4e-01
   5.8e-01   8.4e-01
   4.9e-01   2.7e-01
   7.0e-02   8.0e-01
   9.3e-01   8.7e-01
   4.4e-01   8.6e-01
   3.3e-01   4.2e-01
   8.9e-01   9.0e-01
   4.9e-01   7.0e-02
   9.5e-01   3.3e-01
   6.6e-01   2.6e-01
   9.5e-01   7.3e-01
   4.2e-01   9.1e-01
   6.8e-01   2.0e-01
   5.2e-01   6.2e-01
   7.7e-01   6.3e-01
   2.0e-02   2.9e-01
   9.8e-01   2.0e-02
   5.0e-02   7.9e-01
   7.9e-01   1.9e-01
   6.2e-01   6.0e-02
   2.8e-01   8.7e-01
   6.9e-01   1.0e-01
   6.9e-01   3.7e-01
   0.0e+00   7.2e-01
   8.7e-01   1.7e-01
   6.3e-01   4.0e-02
   3.2e-01   7.3e-01
   4.0e-02   4.6e-01
   3.6e-01   9.5e-01
   8.2e-01   6.7e-01 ];

<span class="comment">% objective point</span>
obj=[0.5,0.5];

figure(1);
<span class="comment">% display the utility function's level sets on some data points.</span>

plot(data(:,1),data(:,2),<span class="string">'o'</span>);
hold <span class="string">on</span>;

[X,Y] = meshgrid(0:.01:1,0:.01:1);
Z=(1.1*X.^(1/2)+0.8*Y.^(1/2))/1.9;

[C,h] = contour(X,Y,Z,[.1,.2,.3,.4,.5,.6,.7,.8,.9],<span class="string">'--'</span>);
clear <span class="string">X</span> <span class="string">Y</span> <span class="string">Z</span> <span class="string">C</span>
hold <span class="string">off</span>;
xlabel(<span class="string">'x_1'</span>);
ylabel(<span class="string">'x_2'</span>);
hold <span class="string">off</span>

m = size(data,1);  <span class="comment">% number of baskets, including 0,1</span>

<span class="comment">% add preference data</span>
Pweak = zeros(m+1,m+1);
<span class="keyword">for</span> i=1:m,
   <span class="keyword">for</span> j=1:m
      <span class="keyword">if</span> (i~=j) &amp; (1.1*data(i,1).^(1/2)+0.8*data(i,2).^(1/2))/1.9 &gt;= <span class="keyword">...</span>
             (1.1*data(j,1).^(1/2)+0.8*data(j,2).^(1/2))/1.9,
         Pweak(i,j) = 1;
      <span class="keyword">end</span>;
   <span class="keyword">end</span>;
<span class="keyword">end</span>;

<span class="comment">% Find consumer preferences</span>
data = [data; 0.5 0.5];
bounds = zeros(m,2);
<span class="keyword">for</span> k = 1:m
    fprintf(1,<span class="string">'Deciding on bundle %d of %d: '</span>,k,m);

    <span class="comment">% Check for u(k) &gt;= u(0.5,0.5)</span>
    cvx_begin <span class="string">quiet</span>
        variables <span class="string">u(m+1)</span> <span class="string">g_x(m+1)</span> <span class="string">g_y(m+1)</span>
        minimize(u(k)-u(m+1))
        subject <span class="string">to</span>
            g_x &gt;= 0;
            g_y &gt;= 0;
            ones(m+1,1)*u' &lt;= u*ones(1,m+1)+(g_x*ones(1,m+1)).*<span class="keyword">...</span>
              (ones(m+1,1)*data(:,1)'-data(:,1)*ones(1,m+1))+<span class="keyword">...</span>
              (g_y*ones(1,m+1)).*(ones(m+1,1)*data(:,2)'-data(:,2)*ones(1,m+1));
            (u*ones(1,m+1)).*Pweak &gt;= (ones(m+1,1)*u').*Pweak;
    cvx_end
    bounds(k,1) = cvx_optval;
    fprintf( 1,<span class="string">'%g'</span>, round(cvx_optval) );

    <span class="comment">% Check for u(0.5,0.5) &gt;= u(k)</span>
    cvx_begin <span class="string">quiet</span>
        variables <span class="string">u(m+1)</span> <span class="string">g_x(m+1)</span> <span class="string">g_y(m+1)</span>
        maximize(u(k)-u(m+1))
        subject <span class="string">to</span>
            g_x &gt;= 0;
            g_y &gt;= 0;
            ones(m+1,1)*u' &lt;= u*ones(1,m+1) + (g_x*ones(1,m+1)).*<span class="keyword">...</span>
              (ones(m+1,1)*data(:,1)'-data(:,1)*ones(1,m+1))+<span class="keyword">...</span>
              (g_y*ones(1,m+1)).*(ones(m+1,1)*data(:,2)'-data(:,2)*ones(1,m+1));
            (u*ones(1,m+1)).*Pweak &gt;= (ones(m+1,1)*u').*Pweak;
    cvx_end
    bounds(k,2) = cvx_optval;
    fprintf( 1,<span class="string">' %g\n'</span>, round(cvx_optval) );

<span class="keyword">end</span>

figure(2);
hold <span class="string">off</span>

<span class="comment">% plot data pt and contour line through it</span>
val = 1.1*sqrt(0.5)+ 0.8*sqrt(.5);   <span class="comment">% value at center</span>
t = linspace(((val-.8)/1.1)^2, 1, 1000);
y = ( (val - 1.1*(t.^(1/2)))/.8 ).^2;
plot(t,y,<span class="string">'--'</span>, [.5 .5], [0 1], <span class="string">':'</span>, [0 1], [.5 .5], <span class="string">':'</span>);
axis([0 1 0 1]);
hold <span class="string">on</span>

<span class="keyword">for</span> k=1:m
   <span class="keyword">if</span> bounds(k,2) &lt; 1e-5,  <span class="comment">% preferred over (.5,.5)</span>
      dot = plot(data(k,1),data(k,2),<span class="string">'o'</span>);
      <span class="comment">%'MarkerSize',8);</span>
   <span class="keyword">elseif</span> bounds(k,1) &gt; -1e-5,  <span class="comment">% rejected in favor of (.5,.5)</span>
      dot = plot(data(k,1),data(k,2),<span class="string">'o'</span>,<span class="string">'MarkerFaceColor'</span>,[0 0 0]);
   <span class="keyword">else</span> <span class="comment">% no conclusion</span>
      dot = plot(data(k,1),data(k,2),<span class="string">'square'</span>, <span class="string">'LineWidth'</span>,1.0,<span class="keyword">...</span>
           <span class="string">'MarkerSize'</span>,10);
   <span class="keyword">end</span>;
<span class="keyword">end</span>;
xlabel(<span class="string">'x_1'</span>);  ylabel(<span class="string">'x_2'</span>);
</pre>
<a id="output"></a>
<pre class="codeoutput">
Deciding on bundle 1 of 40: 0 Inf
Deciding on bundle 2 of 40: -Inf -0
Deciding on bundle 3 of 40: -Inf -0
Deciding on bundle 4 of 40: -Inf -0
Deciding on bundle 5 of 40: -Inf -0
Deciding on bundle 6 of 40: -Inf 0
Deciding on bundle 7 of 40: -Inf -0
Deciding on bundle 8 of 40: 0 Inf
Deciding on bundle 9 of 40: 0 Inf
Deciding on bundle 10 of 40: -Inf -0
Deciding on bundle 11 of 40: 0 Inf
Deciding on bundle 12 of 40: -Inf -0
Deciding on bundle 13 of 40: -Inf -0
Deciding on bundle 14 of 40: 0 Inf
Deciding on bundle 15 of 40: 0 Inf
Deciding on bundle 16 of 40: -Inf -0
Deciding on bundle 17 of 40: 0 Inf
Deciding on bundle 18 of 40: -Inf -0
Deciding on bundle 19 of 40: 0 Inf
Deciding on bundle 20 of 40: -Inf -0
Deciding on bundle 21 of 40: 0 Inf
Deciding on bundle 22 of 40: 0 Inf
Deciding on bundle 23 of 40: -Inf -0
Deciding on bundle 24 of 40: 0 Inf
Deciding on bundle 25 of 40: 0 Inf
Deciding on bundle 26 of 40: -Inf -0
Deciding on bundle 27 of 40: -Inf -0
Deciding on bundle 28 of 40: -Inf -0
Deciding on bundle 29 of 40: -Inf Inf
Deciding on bundle 30 of 40: -Inf -0
Deciding on bundle 31 of 40: -Inf Inf
Deciding on bundle 32 of 40: -Inf -0
Deciding on bundle 33 of 40: -Inf Inf
Deciding on bundle 34 of 40: -Inf -0
Deciding on bundle 35 of 40: -Inf Inf
Deciding on bundle 36 of 40: -Inf -0
Deciding on bundle 37 of 40: -Inf Inf
Deciding on bundle 38 of 40: -Inf -0
Deciding on bundle 39 of 40: 0 Inf
Deciding on bundle 40 of 40: 0 Inf
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="preference_regions__01.png" alt=""> <img src="preference_regions__02.png" alt=""> 
</div>
</div>
</body>
</html>